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MODELS OF BALLISTIC PROPAGATION OF HEAT AT LOW 

TEMPERATURES 

R. KOVACS123 and P. VAN^^s 


Abstract. Heat conduction at low temperatures shows several effects that 
cannot be described by the Fourier law. In this paper, the performance of 
various theories is compared in case of wave-like and ballistic propagation of 
heat pulses in NaF. 


1. Introduction 

The detection and modeling of non-Fourier effects can be vital in cases of energy 
management of very fast processes, in micro-electro-mechanical systems, and also in 
materials with heterogeneous substructures [T]. The second sound - the wave-like 
propagation of heat - is observed at low temperatures in various materials, such as 
NaF, Bi and superfluid He. The phenomenon was predicted by Tisza and Landau 
[HIS] and first measured by Peshkov [3] . Its simplest theory is a memory extension 
of the Fourier equation, the so-called Maxwell-Cattaneo-Vernotte (MCV) equation 
[smE]. Later, the microscopic model of Guyer and Krumhansl [S] proved extremely 
helpful in experimental observations of the second sound in solids, and plays an 
increasingly important role in recent studies [gEoj. The ballistic-type propagation 
is more complicated, to observe, the presence of such effect is very sensitive to 
boundary conditions. Several experiments have detected the latter phenomenon 
[53] , and various theories have been developed to provide explanation and modeling 
[Il|5nilIIl|5g|51|. A short review of the interaction of experiments and theory of 
heat conduction beyond the Fourier law is given in m- 

In this paper, we analyze the low-temperature NaF experiments of Jackson, 
Walker and McNelly, and compare the performance of theories in reproducing the 
experimental data. Our primary concern and reference point is the identification 
of material parameters in our continuum theory, developed in [18] . This theoretical 
framework covers all the known model candidates, including the nine-field equations 
of Extended Thermodynamics, and is based only on universal principles. In |I8j . 
we have used artificial parameters in order to demonstrate the phenomena, while 
in this paper we provide a detailed analysis and modeling. 

2. NaF experiments 

In the experiments analysed here, a heat pulse is applied at the front end of 
the sample, and the relevant measurement data is the rear-side temperature as 
the function of time. In figure one can see several data series for two different 
samples and initial temperatures [5^ . Our focus will be on the middle curve 

on the left-hand side, which shows the characteristic behaviour the best. 
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Figure 1. Original results of NaF experiments: the right one was 
published in [22] and the left one in |23|. L and T denote the 
longitudinal and transversal peaks of ballistic propagation, respec¬ 
tively. 


In what follows, we summarize the available information about the parameters 
and subsidiary conditions of the experiments, which are crucial for theoretical re¬ 
production. 

2.1. Identification of the samples. The papers mention at least six different 
samples, with different sizes and purities. Five of them are mentioned in [23] and 
one in |14| . The crystals can be identihed only by their length and by impurity 
classes such as “very pure”, “less pure”, etc. This twofold identification is crucial 
when we are seeking for the proper parameters of the respective measurement data. 

2.2. The material parameters. Thermal conductivity depends on temperature 
and on the crystal impurity. Walker m presents the thermal conductivity curves 
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for several crystals, the relevant ones being NaF and #2. The peak value of the 
conductivity of is approximately the double of the peak value of the second one. 
This sample has not been used in further measurements, in spite of the seemingly 
better properties. 

Thermal conductivity is not the only relevant material parameter, as one needs 
to use density and specific heat as well. The temperature dependency of thermal 
conductivity is well-discussed in the literature. Bargmann and Steinmann m give 
the value of density and of specific heat according to the Gmelin Handbook [53], 
but those values are different of any of the above-discussed NaF crystals. Hardy 
and Jaswal m have calculated the temperature dependency of specific heat, which 
changes significantly between 9 K and 17 K. This property is analyzed by Coleman 
and Newmann m, too. However, one cannot find information regarding the crystal 
impurity dependency of the specific heat. 

The relaxation times are determined from the measurement via fitting to the 
model equations, hence, different models can lead to different relaxation times. 
Ballistic speed is not an independent further parameter: it can be computed from 
the respective theory once the value of relaxation time is known. 


2.3. Boundary conditions. There are two important questions about the bound¬ 
ary conditions: How long is the input pulse? Can the rear end be adiabatic? 

According to McNelly et al. m the pulse length is between O.lps and Ips. The 
theoretical analyses use various values within this interval. Dreyer and Struchtrup 
apply the lower limit, 0.1/rs m- Yanbao Ma used two different values, 0.24 ps and 
0.36 ps, for the evaluation of the same experiment m- 

Regarding the other question, Jackson et al. |23| do not provide any direct 
information about the rear boundaries. The flattening of the measured curves in 
Figure indicates that the rear side cools back non-negligibly during th eitme scale 
of the measurement. Theoretical modeling did not use this observation. 


2.4. Transversal or longitudinal? The presence of shear and longitudinal bal¬ 
listic propagation indicates that the experiment cannot be modeled in one spatial 
dimension. Moreover, the sketch of McNelly regarding experimental arrangement 
in [28] confirms that, as the specimen is cubelike and the heat excitation affects 
only the middle part of the front face. However, model calculations are mostly 
one dimensional [IIIIIH], except for m- Moreover, usually the transversal, rather 
then the longitudinal, ballistic peak is reproduced in these one dimensional models. 
More properly, something between the two is provided, corresponding to the prop¬ 
agation speed calculated with the one dimensional elasticity parameter, the Young 
modulus. 

On the one hand, it is the longitudinal ballistic propagation that can be con¬ 
nected to one-dimensional models, in the form of volumetric thermoelastic effects 
[36] , while the transversal one cannot be interpreted within such framework. On the 
other hand, the signal speed may be calibrated by the proper material parameters 
in both cases. Hence, in this work we choose the simultaneous reproduction of the 
peak of the transversal ballistic propagation, and of the with the second sound, in 
order to keep our results comparable to results of other theories. 
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3. Non-equilibrium thermodynamics 

In the present study the NaF experiments are modeled in the continuum frame¬ 
work of non-equilibrium thermodynamics. This is the starting point regarding the 
comparison with other models. Therefore, the derivation of our model is briefly 
discussed. Further details are given in [18]. 

Let us consider a rigid conductor, where mass density is constant and material 
time derivatives are identical to partial time derivatives. We start with the balance 
of internal energy, 

ate + V-q = 0, (1) 

where dt is the partial time derivative, e denotes the density of internal energy, V • 
stands for the divergence, and q is the heat flux. The second law reads as 

-b V • J > 0, (2) 

where s denotes entropy density and J is the entropy current density. The as¬ 
sumption of local equilibrium of thermal interaction leads to Fourier’s equation. 
Non-equilibrium thermodynamics can characterize deviation from local equilibrium 
both in the entropy density and in the entropy flux by introducing internal variables 
(dynamical degrees of freedom |32] 1 and current multipliers |17j . respectively. In 
case of heat conduction, it is customary to choose heat flux, q, as internal variable 
[T9l|9]. We also introduce a second order symmetric tensor, Q, as the next approx¬ 
imation. Thermodynamic stability is preserved by assuming a quadratic deviation 
in the entropy density |5^ : 

s(e,q, Q) = Seq(e) - ^q-q-: Q, (3) 

where mi and m 2 are positive material parameters. This choice of field variables 
allows comparison with Extended Thermodynamics. The entropy current is con¬ 
sidered as 

J = b q-fB:Q, (4) 

where b is a 2"'^ order tensor and B is a 3’'*^ order tensor. These quantities are 
called current multipliers - they were first introduced by Nyfri m - and they are 
constitutive functions. Then we can calculate the entropy production: 

9ts-l-V-J = ^b-;^I^:Vq-t-(V-b-mi9tq)-q-|- 

(V-B-maatQ) : Q-bB:VQ > 0. (5) 

Here, I stands for the identity tensor and the triple dot denotes the full contraction 
between the third order tensors. One can separate the thermodynamic forces and 
fluxes according to Table I, and introduce linear relationship between them. 



Thermal 

Extended thermal 

Internal 

Extended internal 

Fluxes 

V • b — mi^tq 

1—1 

1 

V • B — m29tQ 

B 

Forces 

q 

Vq 

Q 

VQ 


Table 1. Thermodynamic fluxes and forces 
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Let US consider an isotropic material and one dimensional propagation. Then 


the linear solution of inequality (§ is 

midtq - d^b = 
m2dtQ - dxB = 



B = 


-kq, 

(6) 

-kiQ -1- ki2dxq, 

(7) 

—k2iQ -f k2dxq, 

(8) 

ndxQ • 5 

(9) 


where the following conditions are required on the coefficients k 2 , ki 2 , n: 

h > 0, ki > 0, k 2 > 0, n> 0, K = kik 2 — ^ 12^21 > 0. (10) 


Eliminating the current multipliers b and B by substitution, evolution equations 
for q and Q are obtained. Then we define the relaxation times Tq = tq = 
the Fourier heat conduction coefficient A = and material parameters K 21 = 

ki 2 = ^ . In order to obtain a system that is close to the hyperbolic one 
of Rational Extended Thermodynamics we simplify the system introducing 
the condition n = ^2 = 0. The resulting system of equations has been called 
ballistic-conductive model [18j . 

The front side boundary condition is a heat pulse modelled by the following 
form: 


qo{t) = q{x = 0,t) = 


Qma 


^1 — cos 

0 


2^-f 


if 0 < t < tp, 


if t > t 


p- 


At the rear side of the specimen, x = L, we consider adiabatic insulation, there¬ 
fore, q{x = L,t) = 0. Initially the fields are homogeneous, and the initial conditions 
are T{x, t = 0) = Tq and q{x, t = 0) = 0. 

After these preparations we introduce the dimensionless variables t, x, T, q, Q for 
time, space, temperature, heat flux and the internal variables, respectively. 


« at 

with 

f T-To 

Tfinal — To 

with 

q 

q = — 
qo 

with 

Q = xl^qoQ- 
V ^21 



y 

pc’ 


X = 


X 

L’ 


rp 7^ I . 

^ final — J 0 H-F ) 

pcL 


1 

90 = — / qoit)dt] 
tp Jo 


The following dimensionless parameters are introduced, corespondingly 


ta = 


(y.tp 

u 


OiTq _ ^ _ ^^12^12 

U’ Z 


( 11 ) 


( 12 ) 


Finally, we obtain the ballistic-conductive equations in a nondimensional form 


fAd{r + d,iq = 
Tqd^q -f <7 -f ta^sT -f kds;Q = 
TQdfQ kd^q = 


0 , 
0 , 

0 . 


( 13 ) 
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Here, all parameters are nonnegative, according to the inequalities (10). This set 
of equations is similar to the set of equations that one can derive from the kinetic 
theory of phonons |T5] . 

It is instructive to eliminate q and Q to reach the partial differential equation of 
temperature: 


+ TQ)diif + d(r - (fg + = 0 


(14) 


One can see that the boundary conditions for this single partial differential equation 
are not straightforward. 

The phase speed is calculated with the help of dispersion relations or hyperbol- 
icity conditions. For instance, the ballistic propagation speed, c, reads 


c = 



(15) 


In what follows, we will use this form to adjust the appropriate propagation speed. 


4. Theory versus experiments 

In this section, we reproduce the results of NaF measurements of Jackson and 
Walker at 13 K (Figure left, middle curve) in the framework of the previously 
outlined continuum theory of ballistic-conductive system of equations, (13), and 
compare the parameters and modeling results with the performance of the 9-field 
approximation of the kinetic theory by Dreyer and Struchtrup, and and the complex 
viscosity calculations of Yanbao Ma [21]. Our calculations have been performed 
with the finite difference method described in 1181. 


4.1. Results of Dreyer and Struchtrup. The field equations of Extended Ratio¬ 
nal Thermodynamics are originated in kinetic theory and modelled through phonon- 
phonon interactions usi dUEsmi]. The 9-held approximation results in a set of 
equations that is similar to ours. The parameters are specific heat, heat conductiv¬ 
ity, and only one of the relaxation times was fitted to the measurements. However, 
their goal was only a qualitative reproduction of the measurement data: the ballis¬ 
tic propagation speed could not match because of theoretical reasons. In (13), we 
have one more material parameter, it, allowing a more accurate modeling of wave 
phenomena. 

Let us consider Figure 5.1/b from [TT]. One needs to consider the material 
parameters of the pure crystal at 13 K to reproduce this simulation. Let the thermal 
conductivity be A: = 13500 based on Figure 3/curve B of |2S], which corresponds 
to the ’’less pure” case. Density is p = 2866 ^ and specific heat is c„ = 1.8 
[ID, where we considered the temperature dependency from |30j . From these, 
thermal diffusivity can be calculated: a = = 2.617^. In [TT], there is no 

explicit data about the length of the specimen, however, there is a correlation on 
the top of the figures: ^ = 23, where x is the sample length, c denotes the Debye 

speed and At is the pulse length. Here, At = 10“^ s is given, c = 


one can 


find c = 2747.5 because = 10.4At is also given. Then one obtains that the 
length of the sample is a: = 6.38mm. The second relaxation time is also given as 
T = 2.lAt. Hence all parameters are reproduced to calculate the dimensionless 
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Rear side temperature 



Figure 2. NaF experiment: reproduction of the Dreyer- 
Struchtrup model, Figure 5.1/b in [TT]. 


coefficients of our model, which turn out to be: 

Tq = 0.0686; TQ = 0.0138; ta = 0.0066; k = 0.102. 

Applying adiabatic boundary condition at the rear side of the sample calculation 
of the time dependency of the rear side temperature results in Figure 

4.1.1. Discussion. Let us note that the length of the crystal does not match with 
the data from [^, where L = 7.9mm. Moreover, the pulse length, which is not 
given exactly in the experimental papers, is different from the one used by Yanbao 
Ma [16| . In their calculations, Dreyer and Struchtrup used Laplace transformations, 
with environmental temperature at the infinity m- This assumption is a particular 
realization of cooling at the rear side of the finite sample. In case of adiabatic 
conditions, the final temperature depends on the energy content of the pulse and 
corresponds to T = 1. Therefore, as a summary, we can say that the difference in 
the material parameters of the two models (the extra k coefficient in our case) does 
not seems to play a role in the qualitative reproduction of their calculation. 

4.2. Results of Yanbao Ma [161 [24]. Yanbao Ma applied the hydrodynamic the¬ 
ory, where relaxation effects are modeled by complex viscosity |35l|T2|. He obtained 
a good quantitative matching with the experimental results including both ballistic 
speeds m- In his case, the duration of the input heat pulse is 0.24ps, and the length 
of the sample is L = 7.9mm. He does not provide information about the rear side 
boundary condition. The published relaxation times are — 0.937 ps; t = 0.248 ps 
and the speed of sound is also given, c = 3186^, which is determined from the 
original measurement data. These papers do not contain information about the 

other material parameters, thus we use the same values as in the previous case, 
2 

a = 2.617 Therefore our dimensionless parameters are 

f, = 0.0393; fQ= 0.0091; ta = 0.01; «: = 0.123. 
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The calculation results in the time dependence of the rear side temperature as 
shown in Figure 


Rear side temperature 



Figure 3. NaF experiment: reproduction of the Yanbao Ma 
model in m- 


4.2.1. Discussion. Apparently, these parameters do not reproduce the measure¬ 
ment data; the difference is significant. The structure of the equations is not the 
same, and the uncertainty in the boundary conditions also seems to be an important 
factor. 

4.3. Results of nonequilibrium thermodynamics m- 

4.3.1. Modification of the model. According to the inherently two-dimensional ar¬ 
rangement of the measurements with a point-like heat pulse at the front end |28j . 
it is reasonable to assume that the conservation of the internal energy is violated. 
Moreover, even if we assume that the thermal insulation at the rear end is not 
perfect, we cannot obtain sufficient heat loss to explain the diffusive tail of the 
measurement in [20) . Therefore, we introduce a source term, a heat exchange be¬ 
tween the bulk crystal and the pulse channel in the balance equation Q in the 
following form: 

dte + d^q = -a{T-To), (16) 

where a is a heat exchange parameter, and Tq = 13K is the environmental tempera¬ 
ture. This is an emulation of the real three-dimensional effect in a one-dimensional 
framework. 

The dimensionless form of the equation is 
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4.3.2. Material parameters. Comparing the information from |28j and [23], the heat 
conduction values seem to be contradictory. The crystal in the experiment accord¬ 
ing to [53] is the one labelled as ”B” in Figure 3 of [53], which has the peak thermal 
conductivity around 15000 However, McNelly points out regarding Figure 32 
in [28] that this crystal is a different one (namely #7204205W [28]), with peak 
thermal conductivity value around 12000 Thus we apply here A = 10200 
instead of the value 13500 based on the Appendix of thermal conductivity 
data in |28j . The relaxation times were determined from the transversal ballistic 
propagation speed and by the proper fit to the second sound peak, and they are 
Tq = 0.355 ps and tq = 0.21 ps. Our additional parameters are n = 0.972mm and 
the heat exchange coefficient is a = 3.2 . 

The dimensionless parameters are, accordingly 

f, = 0.0393; tq= 0.0091; ta = 0.01; k = 0.123; a = 0.149. 

The corresponding solution of the partial differential equation is shown in Figure 
[^ displayed together with is also combined with the originally measured curve 
(Figure [5|. 


Rear side temperature 



Figure 4. NaF experiment: reproduction of the NaF experiment 
at 13K in [20]. The ratio of amplitudes is calculated without con¬ 
sidering the longitudinal wave. 

Due to the two additional parameters, the slope of the diffusive tail and the ratio 
of the amplitudes can be calibrated, too. 

5. Summary and conclusions 

The ballistic-wave-like propagation measurements in NaF by Jackson, Walker 
and McNelly are considered as the best tests of non-Fourier heat conduction at low 
temperature. There are several calculations that qualitatively reproduce the time 
dependence of the rear side temperature. Quantitative reproduction is not an easy 
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ARRIVAL TIME (^sec) 

Figure 5. NaF experiment: combination of the original figure pO] 
and the simulation (Figure]^ . 


task because of the uncertainties of the material parameters and of the boundary 
conditions. 

In this paper, we have compared an experiment at 13K and three theories in the 
framework of our continuum model, published in |18j . We found that this theory can 
be fitted well against the original measurement data, however, it requires different 
parameters than the theory of Extended Rational Thermodynamics or complex 
viscosity hydrodynamics suggest. Moreover, the quantitative reproduction of the 
diffusive tail, considering also realistic boundary conditions at the rear end, requires 
an additional parameter that characterizes internal cooling. This effect seems to 
provide a realistic emulation of the real, three-dimensional phenomena with a point¬ 
like thermal excitation |28j . However, without the anisotropic heat conduction 
properties of the crystals, it is impossible to obtain a complete validation. 
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